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Abetraot. 

This paper describes a solution to the flexural analysis 
of a beam by digital oomputer. Subroutines were written 
which calculate shear, moment, slope and deflection at any 
specified point on a beam if the loading is known. Shear 
defleotion may be included for many support conditions. 
Neither elastic supports nor column effeots are considered. 
The subroutines solve problems of variable beam oross-seotion 
and loading by utilizing an extension of "McCaulay's Method", 
a generalized step funotion. 

The main program provides input-output facilities and 
also solves for indeterminate reaotions on the beam. 

Problems and their solutions, showing the versatility 
of this program, are also inoluded. 
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1. Introduction. 

The routine solution to engineering problems by digital 
computer appears to be limited only by the scope of the 
available programs. A flexible program permitting the solu- 
tion of routine problems in flexural analysis of beams would 
seemingly fill a need in the field of digital computer ap- 
plications to Mechanics. The present work is devoted to suoh 
an idea. 

Five subroutines, "LOAD", "SHEAR", "MOMENT", "SLOPE", 
and "DEFLECT", were written which will calculate the load, 
shear, moment, slope, and deflection, respectively, at any 
indicated point on a beam if the indeterminate quantities 
are calculated elsewhere. The use of generalized step functions 
in this project allows for piecewise linear variation of 
distributed loading, bending compliance, and shear com- 

pliance, 

The main program, "BEAM3", first generates and then 
solves a set of simultaneous linear equations, calculating 
the indeterminate quantities of the beam. The above mention- 
ed subroutines are then utilized to calculate the remaining 
flexural quantities. This program also provides input-output 
facility. 
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2. Definition of the problem. 

The problem of structural beam analysis is one in 
which the external loading applied to the beam is well de- 
fined, as is the geometric configuration and the nature of 
the beam itself. Where information about such a beam is re- 
quired, this information will most likely fall into the 
following categories: 

1. Internal shear forces, V, and moments, M. 

2. Slope,©. 

3. Deflection, y. 

4. Redundant reactions, either forces, Ry, or moments, 

RM Z . 

"Flexural analysis" is used here as meaning the process by 
which this information is generated from the known data. 

In attacking this problem of flexural analysis, the 
intention was to limit the scope as little as possible. 
However, several customary simplifications and limitations 
were made. First, elastic supports and beam column effects 
were not considered. Limitations concerning the type and 
number of external loads and the geometry of the beam will 
be discussed later. Also the basic relations which follow, 
inherently restrict the deflection of the beam to small values 
and the internal stresses to values below the elastic limit 
of the structural material. Shear deflection, which is 
generally neglected in a problem of this type, has been 
Included in some problems, 
a. Bending deflection (y^) . 

The Bernoulli-Euler-Navier equation governing the elastic 
curve of a beam due to bending, consistent with the sign 
convention shown in Appendix B is: 



2 



« curvature * (d 2 y/dx 2 ) [l+(dy/dx) 2 ] ^ 

2 2 

and we use the approximation, curvature d y/dx . 

The following equations, derived from the above equation 
and shown in Timoshenko /5/,^ with differences due to the chosen 
sign convention, describe the flexure of a beam (See Appendix B 
for notation). 



W * q(x) , 
V * f\ Wdx , 
M * _/Vdx , 




yi - J r 9 1 dx , 
b. Shear deflection (y 2 ). 

The equation governing the shear deflection of a beam, 

also given in Timoshenko /5/» compatible with the stated 

sign convention becomes: 

dy~ KV 

dx AG 

It is evident that: 




o. Total deflection (y). 

The equations used in this thesis are the result of the 
sum of the above shear and bending deflections. They are: 

W = q(x) , 



Numbers so indicated refer to the bibliography. 
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Note that step-discontinuities may appear in all of the 
above equations, except in the equation for deflection. These 
discontinuities limit the methods that can be used in handling 
the equations. 

Three methods that could be used to solve the equations 
are a pure analytical approach, a step-by-step numerical 
integration scheme, and the "curly bracket" method. Of 
these methods, the analytical approach is generally restrict- 
ed to a specific problem and is unsuited for adaptation to 
digital computer programming. The numerical integration scheme 
can be used in the computer solution of flexural analysis 
and may have been programmed; however, our search of literature 
has not uncovered any such program. The "curly bracket" 
method was chosen for use in this thesis because it was felt 
that it would result in a more exact solution requiring less 

*In these equations, 2^ and£M represent the sum of 
the point forces and moments, respectively, applied to the 
beam. 

**C^ and C 2 are constants of Integration. 
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computer time, particularly if results were required at only 
a few points. 
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3. "Curly bracket" method 

The "curly bracket" method is based upon the properties 
of a convenient notation defined as follows: 

{x-a} n = 0 i f x <a 

= ( x- a ) n if x> a 

where n is an integer^O. 

Appendix A contains further properties of expressions contain- 
ing curly brackets. The notation is generally credited to 
McCaulay /l/, but similar, though less useful notations can 
be traced to A.Fdppl and St. Venant. Also, the ideas in- 
volved are related to the use of the Laplace Transformation 
with the result that this method or approximations to it are 
reinvented frequently /&/ /9/. 

Using this method allows one equation, valid over the 
entire length of the beam, to be expressed for each of the 
basic equations of flexure. The particular value of this 
method is that it permits the integration of these equations 
to be carried out easily and has the virtue of adjusting the 
constants of integration most conveniently. 

In addition, for use with a digital computer, the "switch- 
ing" property of the curly bracket can be easily handled in 
FORTRAN by use of an "IF" statement. 

A simple example should convince most readers of the 
convenience of this method. Examine the following problem, 
which consists of a uniform beam, simply supported at each 
end with a concentrated load at the center of the span. 
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Fig ♦ 1 

['he following equations lead to the “ solution of this 



oroblem. 
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The only unknown value, Cn, in the equations can easily be 
calculated, using the boundary condition, y = 0 at x-2a. 



C, 



-8E 



Thus the resulting equation defining deflection at any point 



uor.i 0- 



y 



x 3 

, t?i_ 
" 1 9 



e. .--.3 

v t :""- 1 



As hcC-ulay used this method, it was suitable only for 
beams of uniform cross-section. Quinlan /3/ suggested - way 
of expanding the net hoc. so include beams with piecewise ccn- 
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stant variation of bending compliance, In this thesis 

the method has been further expanded to include beams with 
piecewise linear variation of both bending compliance and 

, K 

shear compliance, . 

AG 
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4. Implementation. 

Although the formulas shown in Appendix A could be ex- 
panded further, a piecewise linear variation of distributed 
loads, bending compliance, and shear compliance was consider- 
ed capable of providing sufficient accuracy and flexibility 
for the resulting program. 

Fig. 2(a) Illustrates a distributed load on a beam 
whose origin is at x=0. Fig. 2(b) shows an element of the 
distributed load. 





Fig. 2(b) 



The analytical expression for this element is: 



Q(x) - QaJx-A} 0 + ^^(x-A} 1 - ^{x-B] 0 - {x-B> 1 



The same relationship Is used for an element, f(x), of bend- 
ing compliance and an element, g(x), of shear compliance upon 
substituting f for Q, and g for Q, respectively. 

The basic equations of flexure can now be expressed In 
the same manner as they were utilized In programming. They 
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are; 



NOL 

W(x) - Z Q,(x) 
1=1 1 



(1) 



£ NOF . NOR 

W(x)dx + £ F y fx-x f ] + Z R „ f x “ x r 1 

i«l y i L i - 5 i*l y l* i J 



( 2 ) 



where F are upward, known forces at x^x*. and R are up- 
y l 1 y i 

ward, unknown reactions at x^x 



NOM 



NORM 



M(x) = \ V(x)dx + Z M, f x “ x m 1° + Z R M £x-x V 

Jo i=l z i^ V i-1 z i* rm ii 

(3) 

where M are clockwise, known moments at x-x and RM„ are 
z i m i z i 

clockwise, unknown reactive moments at x«x 



rnu 



G 



dy = 
dx 



NOI 

M(x) Z f<(x)dx - 
i=l 



NOK 

v(x) Z Si( x ) 

i=l 1 . 



♦ c- 



(4) 



where is the constant of integration. 

Appendix A indicates how the product and the integration of 
the product of the curly bracket symbols are formed. 



^x r x r _ NOI 

M( 




1NUJ. "I fX~ NOK "I 

1)2 f i( % )J dx v (x) Zgj(x)dx + + C 2 

i - 1 - 1-1 ( 5 ) 

where C 2 is the constant of integration. 

Equations (1) through (5) when properly coded in FORTRAN 



The numbers, NOL, NOF, NOR, NORM, NOI and NOK, appear- 
ing above the summation signs are program input parameters 
indicating the number of terms for each problem. 



10 



form the core of the five basic subroutines? LOAD, SHEAR 
MOMENT, SLOPE and DEFLECT. It must be noted before discuss- 
ing the subroutines that equations (2) through (5) Involve 
the unknown values C 2 , R y , and RM z , where there may be 
as many as ten values each of R y and RM Z . The calculation of 
these values is done in the main program BEAM 3 and will be 
discussed in section 6. 
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5. The basic subroutines 



Subroutines LOAD, SHEAR and MOMENT, 

The above named subroutines, as their names imply, are 
employed to calculate the total distributed load, shear and 
moment, respectively, at any desired point on the beam. The 
expanded forms of equations (1), (2) and (3) were used in 
writing these subroutines. Utilizing the distance, x^, as 
the entering argument, the subroutines select and sum the 
appropriate terms in the corresponding equation by inspect- 
ing all loading terms and discarding those in which the value 
within the curly bracket is negative. Actually the load, 
shear and moment are calculated a small increment, £ - 10 , 

to the right of x^ in order to obtain the right side value 
of any step-discontinuities occurlng at the point, x ^ (See 
flow charts on pages 13 , 14 and 15 ). 

Subroutine SLOPE. 

This subroutine is equation (4), excluding the constant 
of integration, C^, in FORTRAN language (See flow chart on 
page 17). Note that the equation can be expressed as a 
summation of terms, which are constants times the multiplica- 
tion of two step- functions of the form C £x-aj n £x-b} m , or the 
summation of the integral of the same type of terms, where n 
equals either 0 or 1. 

Function subroutines, FTON and FT1N (see Appendix A), 
will calculate the value of £x-a} ° [x-bj m and [x-a^ ^{x-b} n , 
respectively. Function subroutines, ENTON and ENT1N (See 
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Flow Chart of LOAD 




Note: Numbers preceded by a number sign, refer to 

the statement number in the listing of LOAD. 
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Flow Chart of SHEAR 




Note: Numbers oreceded by a number sign, refer to 

the statement number In the listing of SHEAR. 
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Flow Chart of MOMENT 








RETURN 



Note: Numbers preceded by a number atgr, #, refer 

statement number In the listing of MOMENT, 
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Appendix A), will calculate the first integral of these 
functions. Equation (4) can now be expressed as the summation 
of constants multiplied by the appropriate function subroutines. 

SLOPE calculates all combinations of M(x) and f(x) by 
stepping through all the f(x) terms and varying all the terms 
in the moment equation at each step. It then steps through 
all the g(x) terms varying all the terms in the equation for 
shear at each step. The sum of these terms is the value of 
the slope, minus the integration constant, 0^, at the entered 
distance, x. 

Subroutine DEFLECT. 

At this point the author wishes to point out the similar- 
ity of the equations for slope and deflection. Since the 
integration constants, C-^ and C 2 , are again excluded from the 
subroutine, and the function subroutines, DITON and DITIN 
(See Appendix A), calculate the seoond integral of the mul- 
tiple step- functions, subroutine DEFLECT is identical to 
SLOPE except for substituting ENTON for FTON, ENT1N for FT1N, 
DITON for ENTON and DITIN for ENT1N. The result is a sub- 
routine that will calculate the deflection, minus the con- 
stants of integration, at any required point on the beam. 
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6. Program BEAM3. 

The main program, BEAM3, can be divided into four dis- 
tinct steps. These steps are: reading the input data, cal- 

culation of indeterminate quantities, calculation of the 
flexural quantities, and the output of the solution. 

The program first reads the input data which includes 
the length of the beam, all external loading, the location of 
all reactions, the flexural properties of the beam, and a group 
of parameters giving the number of each type of load or re- 
action. It also provides for reading a predetermined value 
of deflection, PDEFT, at each reaction and a predetermined 
value of slope, PSLP, at each moment reaction. 

The program then proceeds to solve for the indeterminate 
quantities of the beam by utilizing the two equations of 
static equilibrium plus equations obtained from the known 
boundary condition at each reaction or reactive moment. The 
first two equa;ions are simply the summation of the forces 
on the beam equals zero and the summation of the moments 
about the right end of the beam equals zero. The remainder 
of the equations are obtained from the fact that the slope 
at each reactive moment or the deflection at each reaction 
is known. Using this fact in the appropriate one of either 
equation ( 4 ) or ( 5 ) results in a set of linear equations for 
the indeterminate quantities. 

The equations thus obtained are generated in a matrix 
form that is compatible with subroutine GAUSS2 (See Appendix E) . 
This subroutine is then called to solve the equations simulta- 
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neously. If the matrix is singular, GAUSS 2 exits to the call- 
ing program which in turn prints "MATRIX SINGULAR" and then 
stops. If not, the subroutine solves the equations and re- 
turns the calculated values to the main program. The most 
serious limitations of this program are imposed by the gen- 
eration or solution of the simultaneous equations, therefore 
it is appropriate to include and explain several of them at 
this point. 

First, at least one reaction, even if it is later forced 
to equal zero, must be included with the input, or a singular 
matrix will be generated. 

Including shear deflection in the solution must also be 
greatly limited at this point. Since in general the shear 
includes step discontinuities, they will also appear in the 
slope if shear deflection is included. The program generates 
an equation for slope at each moment reaction. If by chance, 
either a reaction or point force occurs at the same point as 
a moment reaction, the slope at this point is double valued. 
The program is unable to distinguish which of the two values 
is required except when it occurs at the extreme ends of the 
beam. For example a beam "built in" at both ends is accept- 
able. If the beam is "built in" at only one end, it must be 
solved in a manner such that the "built in" is at the origin 
of the beam. 

The program then divides the beam into a number of 
increments, NOP-1,* and calls the subroutines, LOAD, SHEAR, 

*NOP is another parameter read in at the start of BEAM 3 , 
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MOMENT, SLOPE and DEFLECT to calculate the flexural quantities 
at each of these increments. 

The output consists of printing the distance from the 
left end of the beam to the point in question and the values 
of the load, shear, moment, slope and deflection at that 
point. The values of all of the indeterminate quantities are 
also printed out. 
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Fig* 7 Simplified Flow Chart of BEAM3 
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7. Testing of Program. 

The testing of the program and the subroutines could, 
in theory, continue indefinitely because of the unlimited 
number of beam configurations. The author has attempted to 
show (See Appendix D) solutions to test problems of simple 
configurations and constant bending and shear compliance to 
compare with classical solutions. Then a statically deter- 
minate problem with a varying bending compliance was solved 
and compared to a solution obtained from a numerical in- 
tegration scheme. The most rigid test for the program was 
in constructing a problem for a beam of unlikely configuration 
(See Problems 6 and 7 of Appendix D), and solving the flexural 
analysis of this beam from both ends. Note that the only 
significant difference between the two solutions is the opposite 
signs for shear and slope, which would of course, be expected. 



22 



8. Conclusions and Recommendations. 

The thesis shows a method of including piecewise linear 
variation of both bending compliance and shear compliance in 
the flexural analysis of beams by extending "McCaulay's Method". 
Adapting this idea to the digital computer resulted in a 
program capable of analyzing a very general type of beam. 
However, had time permitted, several changes would have been 
made to the program. First, internal monitoring of the input 
data would be an asset to any user. Also, using the same 
general method of attack, this program would have been ex- 
panded to include bending compliance and shear compliance of 
a piecewise parabolic nature. 

Other additions to the program, much larger in scope, 
would be including elastic supports and/or beam column effect. 

It would be interesting to compare the results of this 
program to those obtained from a program utilizing the numer- 
ical integration method. Also, a more exhaustive study of the 
program could be made to determine accuracy of, and time 
requirement for, solutions. 
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APPENDIX A 



"Curly bracket" method, definitions and formulas 
We define: 

1. [x-aj n = (x-a) n , x>a 
= 0 , x < a 

where n is an integer^ 0 

From this definition, it is possible to establish the 
following formulas. Although it would be possible to express 
these results in terms of curly brackets, in most cases the 
resulting expressions would be rather complicated. The forms 
exhibited below have the advantage of being most readily 
conformable to the decision commands available in FORTRAN. 




(x-a) 

n+1 



n+1 



0 



x> a 

Otherwise 



3. [x-a] 0 ;x-b3 n 

k. ^x-aj-^x-b} 11 

5. f [x-a) 0 {x-b} n d x 



( x-b ) n 
0 



x> a and x>b 
Otherwise 



(x-b) n+1 + (b-a)(x-b) n , x > a 

and x >b 



(x-b) 

n+1 



A*=*l 

n+1 



n+1 



Otherwise 
n+1 



(a -b) 

n+1 



n+1 



x >a >b 



x>b>a 



0 



, Otherwise 
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6. f £x-a} 1 [x-bj n dx 

)o 



n+2 



(x-b) 
n+2 

, v »\ nf 2 

( a-b ) 

(n+1) (n+2) 

(x-b) + 

n+2 



+ (b-a) (x-b) 

n+1 

, x>a>b 
(b-a) (x-b) n+1 



n+1 



n+1 



x>b>a 



Otherwi se 



rr B- a 3°^- b 3 n ^ dx 

(x-b)"* 2 
(n+1 



n+1) (n+2) 

, , .n+2 

(a-b) 

(n+1) (n+2) 

, v n+2 
(x-b) 

(n+1) (n+2) 
0 



(x-a) (a-b) 
(n+1) 



n+1 

x >a >b 

x >b ;>a 
Otherwise 



»• a £-*3 i £- b r^ d * 



(x-b) 



n+3 



(n+2) (n+3) 



[ »-a, (a-b) 
+ ( 1 ) n+2) 

t i^n+3 
( x-b) 

(n+2) (n+3) 



n+2 



n+2 



(b-a) (x- b) 
(n+1) (n+2) 

2<a-b>"* 3 



(n+1) (n+2) (n+3) , 

n+2 



( b-a) (x- b) 
(n+1) (n+2) 



x>a >b 
x> b >a 



, Otherwise 



Function Subroutines FTON, FT1N, ENTON, ENT1N, DITON, DITIN 
The above named function subroutines were written and 



utilized to evaluate each of the preceeding formulas numbered 
3 through 8. The values x, a, b and n were used as entering 
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arguments to these functions. The function then became the 



variable as 


follows: 


FTON 


= |x-a]Vb} n 


FT IN 


* [x-a} 1 £x-bj n 


ENTON 

ENT 1 N 

DITON 

DITIN 


- C £x-a] 0 [x-b] n dx 

SO 

- f ^ £x-b} n dx 

JO 

- , 00 S' a i° d * 

’ 0l£- a i 1 fc _b ? nd?dx 



As examples, flow charts for FTON and ENTON are shown in 
Figs. 8 and 9. 
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Fig. 8 



Flow Chart of FTON 
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Fig. 9 Flow Chart of ENTON 
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APPENDIX B 



Sign Convention and Notation 

Sign Convention 

The following sign convention has been employed: 

y 

t 

i- x y * deflection 




O = slope = 

dx 



»(r i) 
v tczzu 



M » bending moment 



V = shear force 



W 



TTTTT 



¥ = distributed loads 



Note: 1. Point forces, F y , and point reactions, R , are 

positive in the positive y direction. y 

2, Point moments, M a , and point reactive moments, 
RM 2 , are positive in a clockwise sense. 



Notation: 

THESIS 

W 

V 

M 

e 

y 



PROGRAM 

C 

W 

V 

AM 

Y1 

Y 



Length of beam 

Total distributed load at a point 

Shear 

Moment 

♦ 

Slope 

Deflection 
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F y 


FY 


Point load or force 


R y 


RY 


Point reaction 


M z 


AMZ 


Point moment 


RM Z 


RAMZ 


Point reactive moment 


Q(x) 




Element of distributed load 


Q a 


AQ 


Distributed load (left end) 


% 


QjB 


Distributed load (right end) 


f(x) 


EA 


Element of bending compliance 

(It) 

Bending compliance 






(left end of element) 




EB 


Bending compliance 
fright en<? of element) 


g(x) 




Element of shear compliance 

(K_) 
' AG ; 




GA 


Shear compliance 
(left end of element) 




GB 


Shear compliance 
(right end of element) 


X f 


XFO 


Distance to Fy 


x r 


XR 


Distance to Ry 


Xm 


XM 


Distance to M z 


x m 


RXM 


Distance to RM Z 


A g 


AG 


Distance to G a 


B g 


BG 


Distance to G^ 


A e 


AE 


Distance to E a 


B e 


BE 


Distance to E^ 


A 


A 


Distance to Q. 

cl 


B 


B 


Distance to 


NOL 


NOL 


Number of distributed loads (25) 



3 

Numbers following explanation of terms NOL--NOK indicate 
the arbitarily chosen maximum value of each for this program. 
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NOF 


NOF 


Number of point forces (100) 


NOR 


NOR 


Number of point reactions (10) 


NOM 


NOM 


Number of point morasnts (50) 


NORM 


NORM 


Number of point reactive 
moments (10) 


NOI 


NOI 


Number of bending compliance 
elements (25) 


NOK 


NOK 


Number of shear compliance 
elements (25) 


PDEFT 


PDEFT 


Predetermined deflection 
at a reaotion 


PSLP 


PSLP 


Predetermined slope at a 
reactive moment 


NOP 


NOP 


Number of points at which 
the flexural quantities will 
be calculated. 
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APPENDIX C 



General Information Concerning Use of Program 

The composite program, BEAM 3 plus subroutines, was 
written in FORTRAN-60, compiled and tested on a Control 
Data Corporation l6o4 digital computer. It consists of 412 
FORTRAN statements requiring about 580 cards. The computer 
storage requirement is about 12,500 cells, however, this 
number can be reduced by changing the dimension statements 
which would of course, reduce either the maximum number of 
external loads, elements of shear compliance or bending compli- 
ance, or number of points at which flexural quantities could 
be calculated. 

The FORTRAN statements used in this program are: GO TO, 

computed GO TO, IF, STOP, DO, CONTINUE, FORMAT, READ, PRINT, 

WRITE OUTPUT TAPE, FUNCTION, SUBROUTINE, RETURN, CALL, DIMENSION, 
COMMON, and END. 

There are no conversion constants incorporated in this 
program. As such all input must be in a oonsistant set of 
units. 

"FORMAT" for Input Data 

The input data for a problem is read by the program 
with a FORMAT of 14 for all "fixed point" variables and F20.0 
for all floating joint variables. The arrangement of the 
data on the input oards is shown as follows: 
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Group 

no. 


No. of cards 
in group 


Fields 


Entries 


1 


1 


F20.0 


G 


2 


1 


814 


NOP, NOL, NOF, 


3 


NOL 


4F20.C 


NOM, NORM, NOI 
A, B, A &> A-^ 


4 


NOF 


2F20 . 0 


Xf, Fy 


5 


NOR 


3F20.0 


x r , 0.0, PDEFT 


6 


NOM 


2F20.0 


V M 2 


7 


NORM 


3F20.0 


x m , 0.0, PSLP 


8 


NOI 


4F20 . 0 


A e* ®e» E a» 


9 


NOK 


4F20.0 


A g» B g* & a* °b 



As an example the Input data cards for test problem #1 
In Appendix D was as follows; 



CARD 


Is 


100 






CARD 


2; 


21 1 0 


2 0 


2 10 


CARD 


32 


0 . 


100 


-500. -500. 


CARD 


4s 


0. 


0 . 


.0 


CARD 


5? 


100 


0 . 


oO 


CARD 


6s 


0 . 


0 . 


0 . 


CARD 


72 


100 


0 . 


0 . 


CARD 


8s 


0 . 


100. 


.000000001 .000000001 
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APPENDIX D 



Test problems and solutions 

Test problem 1 

This problem consists of a beam of uniform cross-section, 
"built-in" at both ends, and loaded with a uniformly dis- 
tributed weight of 500#/in. (See Fig. 10) 



X 

X. 
X 
v ~ 

x- 

X 

X" 



1 



Q s 500#/ in. 

LL.1....L 



100 in". 



X 

X 

-x 

x 



1 

El* 

K 

AG 



10 ~ 9 #in. 2 



« 0 



Fig. 10 



Test problem 2 

This problem is identical to problem 1 except that the 
slope at either end of the beam is forced to values other 
than zero. 



0 

left end 



-.002 



p. 

right end 



.002 
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SOLUTION TO TEST PROBLEM 
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SOLUTION TO TEST PROBLEM 
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Test problem 3 

This problem consists of a beam of uniform oross-section, 
"built-in 1 ' at one end and simply supported at the other, and 
loaded with a uniformly distributed weight. (See Fig. 11) 



V 



Q, = lO#/in„ 



znoiErra 



\ 

\ 



V- 



40 in, 



1 

El 

K 

AG 



, 1332 x 10 



-7 1 



#in. 2 



Fig. 11 



Test problem 4 

This problem is identical to problem number 3 except 
shear deflection is included in the solution. 

~~ - .159 x 10“ 6 l/#in. 
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SOLUTION TO TEST PROBLEM 



z Of'OCNiPOcOOCM*— Ojd‘ 0 ' 0'0 0'00 C\iO 

O OCT'OCDOin*— Or'OinOrOrOCOOOunOCXDLO'OO 

— OCOC^O— CMJ* — o*— Or— tnOCOrOOOCMO 

H O^I^OO'rOr-OOOr-rO.=f C 0 r-CT'r-OC 7 'ino 

o oo*— roinoooroinof^coooN-o^fcvjO'Oroo 

OJ OOOOOO^^r-rr-r-r-f-r-r-f-r-OOO 

-J ooooooooooooooooooooo 

CL ooooooooooooooooooooo 

UJ ••••••••««••••••••••• 

o I I I I I I I I I I I I I I I II I I I 



OO^OrOOOO-^rOroO — O' — — OrOrO-^COO 
o oo ro r- c> •— csj tn ^ uo o r- n* •— on o o r>- <5 

uj 0 ' 0 ^^r^c\iaof^oo'^f'Oc\Josjc\iC 7 'Ln< 5 r>opor-. 

a. 0 - 3*000 — CM r ~ 0 C' 0 . 3 ‘»-“ — . 3 ‘f'‘-C N CM- 3 - 0 rwN* 

O OOO— — — — — OOOOOOOO— — — — — 

-J ooooooooooooooooooooo 

00 ooooooooooooooooooooo 

ooooooooooooooooooooo 



I I II I I I I I II I 



ooooooooooooooooooooo 
h- ooooooooooooooooooooo 

2 • • • ♦ • 

ULI ooooo oooooooooooooo 

ST OCMCOCOCM 00 CM CM CO O CO C>JfsJ COO CD CM CM CO 

o oino'oro csiLnr^cooo'-'-oocor-Lncsi 

X CM — I I — »— •— •* — 

I I I 



ooooooooooooooooooooo 
O' ooooooooooooooooooooo 

< 

LU ooooooooooooooooooooo 

X com*— c^r-inro— cr^cnro^^rou^O'— roun 

t/> CMCMCMr— — — | | I I I — — — 

I I I 



o 

< 

o 

-J 



ooooooooooooooooooooo 

ooooooooooooooooooooo 



ooooooooooooooooooooo 



I I I 



II I I I II I II I I 



LU 

o ooooooooooooooooooooo 

2! ooooooooooooooooooooo 

< •••♦♦•••••••••••••••• 

H CSJ^'OCJOOCVJja-'OOOOOJJ-OCOOCsI-^OOOO 

00 — CMCMCMCMCMrOrorni'OrO-d* 



o 

o 

oo • 
ooo 
• «o 

ooo 
cnmcMOo 
CM — I OO 

oo 
oo 
II oo 
—oo 
II II — oo 
— — oo 

— CM*— • • 
Nj 

— -2: ii ii 
>->< — (M 

o£a:a:oo 



a 



39 



SOLUTION TO TEST PROBLEM 4 
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Test problem 5 

This problem consists of a homogeneous circular shaft, 
three inches in diameter at the center and tapered to one 
and a half inches in diameter at both ends. The total weight 
of the shaft is 100# and a uniformly distributed load total- 
ing 200# 1 8 placed on the central section. The shaft is 
supported by a uniformly distributed reaction at each end. 
Points A and B are points of zero deflection. This problem 
was approximated for the program as shown in Fig. 12. 




Fig. 12 



Weight of shaft: 
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SOLUTION TO TEST PROBLEM 
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Teat problem 6 

This hyperstatic problem consists of a beam with a very 
unlikely configuration as shown in Fig. 13 . The bending com- 
pliance varies as shown in Fig. 14. Shear deflection has been 

K 

oml tted, 1 . e . 9 — 0 s 





Fig. 14 



Test problem 7 

This problem is identical to problem number 6 except the 
beam was turned end for end. 
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SOLUTION- TO TEST PROBLEM 6 
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APPENDIX E 



Subroutine GAUSS2 

This routine is a modified form of the GAUSS2 used by 
C. B. Bailey /7/ in his program, LINEQN, for solving simul- 
taneous linear equations by Gaussian elimination. Bailey 
was solving the matrix equation, Ax^B, with a possibility 
of 50 vectors B for each matrix A, but this program requires 
only one vector B with each matrix A. Also since the maximum 
number of equations in BEAM3 is 22, the size of matrix A was 
reduced from 100x101 to 22x23. 

At each step of the elimination, the value of the diagonal 
element is compared to a defined 11 zero ". If it is smaller 
than this quantity, the matrix is considered singular and 
an error output is returned to the calling program. This 
is the "matrix singular" test used by BEAM 3 . 
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Listing of program and subroutines 
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BEAM3 


48 


LOAD 
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MOMENT 
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SLOPE 
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DEFLECT 


61 
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64 


ENTON and ENT1N 


65 


DITON and DITIN 


66 


GAUSS 2 


67 
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PROGRAM BEAM3 

DIMENSION W( 1001 ) , V ( 1 001 ) , AM ( 100 1 ) , Y 1 ( 1 00 1 ) , Y { 1 00 1 ) , XX( 100 1 ) , 

1 A(25) , B( 25) ,QA(25) ,QB(25) ,XFO( 100), FY( 100) ,XR( 1 0),RY{ 10) , 
2XM(50) ,AMZ (50) ,RXMMO) ,RAMZ( 1 0 ) , AE ( 25 ) , BE ( 25 ) , E A { 25 ) , EB ( 2 5 ) , 
3AA(22,23) , XA(22) ,PDEFT(1 0),PSLP(10) 

4 , 1 TITLE ( 10) ,GA (25 ) ,GB(25) , AG (25) ,BG (25 ) 

COMMON A,B ,QA,QB, XFO, FY, XR , RY , XM , AMZ , RXM , RAM Z , AE , BE, EA,E3,N0L, 
1 NOR, NOF,NOM, NORM, NOI , NO K , G A , GB , AG , BG 
READ 70, C 

70 FORMAT (F20.0) 

READ 71 , NOP , NOL , NOF , NOR, NOM , NORM, NOI.NOK 

71 FORMAT (814) 

IF(NOL) 51,51,50 

50 READ 72 , ( A ( I ),B(I),QA(I) , QB ( I ) ,1 = 1, NOL) 

72 FORMAT ( 4F20 . 0 ) 

51 IF( NOF) 53, 53,52 

52 READ 73, ( XFO(N) ,FY(N) ,'4=1, NOF) 

73 FORMAT (2F20.0) 

53 IF( NOR) 55, 55,54 

54 READ 74, (XR(M),RY(M) ,PDEFT(M),M=1 ,NOR) 

74 FORMAT (3F20.0) 

55 I F ( NOM ) 57 ,57 ,56 

56 READ 73, (XM( K) ,AMZ(K) ,K=1 ,NOM) 

57 I F ( NORM ) 59,59,58 

58 READ 74, (RXM(KA) , RAMZ ( KA ) , PSLP( KA) ,KA=1,N0RM) 

59 IF(NOI ) 6 1 ,61,60 

60 READ 72, (AE( IA), BE ( I A) ,E A( IA) ,EB( IA) , I A= 1 , NO I ) 

61 I F ( NOK ) 1 , 1 , 62 

62 READ 72, < AG( IB) ,BG( IB), GA(IB) ,GB( IB), IB=1,N0K) 

1 NORl = NOR+1 

C 1 = 0. 

C2= 0. 

NORM 2 =NOR M + 2 
N0RM3 =NORM + 3 
NOQ = NOR +NORM 
NCQ1 = NOQ + 1 
N0Q2 = NOQ + 2 
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N0Q3 = NOQ + 3 
EPSX= 1 0 . ** ( -8 ) 

DO 2 K= 1 , NCR 

2 AA(1,K) * 1. 

DO 3 K= NCR 1, NOQ 2 

3 AA( 1 ,K ) = C. 

AAA=G. 

IF (NOF) AS ,45,35 
35 DO 4 N = 1 , NC F 

4 AAA = AAA + FY ( N) 

45 I F ( NOL ) 48,48,46 

46 DO 5 1=1, NCL 

5 AAA = AAA + ( ( QA ( I ) +QB ( I 1) /2 . ) * (B(I)-A(I)) 

48 AA ( l,NOQ3)=-AAA 

DO 6 K =1 , NCR 

6 AA ( 2 , K ) = C-XR(K) 

IF{NCRM)9,9,7 

7 DO 8 K=NOR 1 , NOQ 

8 AA { 2 , K ) = 1 . 

9 DO 1C K=NCQ 1 , NOQ2 

10 AA ( 2 ,K ) = 0. 

CALL MOMENT ( C ,A AA ) 

AA ( 2 ,N0Q3 ) =-AAA 
IF ( NORM 117,17,11 

11 DO 16 J=3, N0RM2 
DO 13 K= 1 , NCR 
AA ( J , K ) =0. 

DO 12 I A = 1 »NOI 

AAA = EA( I A ) *ENT 0N( RXM ( J-2 ),AE(IA),XR(K),1) - EB(IA) *ENTON ( RXM ( J-2 ) , 
1 BE( IA) ,XR(K), 1 ) +((EB(IA)-EA(IA))/(Bc( I A )-AE ( I A) ) ) * ( ENT 1 N ( RXM ( J-2 ) 
2, AE ( I A) ,XR(K) , 1 )- ENT 1N( RXM( J-2) , B E { I A ) ,XR(K), 1 ) ) 

12 A A ( J , K ) = A A ( J , K 1 + AAA 

IF (NCK 1 13,13,900 

900 DO 91 C I B= 1 , NOK 

IF ( NORM-1 1902 ,901 ,9C2 
902 IF( J-N0RM21901 ,904,901 
9 C 1 X = RXM ( J-2 ) -EPSX 
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AAA = GA ( IB ) * FTO N ( X , AG ( I B ) , XR ( K ) , 0 ) - GB(IB)* FTOMX 

1 , BG( 18), XR(K) ,0) + { (GB{ IB)-GA(IB) )/{BG(IB)-AG( IB) ) )* ( FT IN ( 

2 X , AG(IB) ,XR(K) ,0)- FT) N ( X , BG ( IB ) , XR( K ) , 0 ) ) 

GO TO 910 

* 

904 X=RXM ( J-2 ) +EPSX 

AAA = G A ( 16)* FTO N ( X , AG ( I B ) , XR { K ) , 0 ) - GB ( I B ) * FTCN(X 

1 , EG C IB) ,XR(K) ,0) + ( (GB( IB)-GA(IB))/(BG(IB)-AG( IB) ))« ( FT 1 N ( 

2 X , AG(IB) ,XR(K) ,0)- FT 1 N { X , BG ( I B ) , XR { K ) , 0 ) ) 

910 A A ( J , K' ) = A A ( J , K ) -AAA 

13 CONTINUE 

IF( N0C-N0R1 ) 1 55, 1 35, 135 
135 DO 15 K=N0R1,NCQ 
AA{J,K)= 0. 

DO 14 I A = 1 ,NOI 

AAA= EA ( IA) *ENTON (RXM( J-2) ,AE( I A ) ,RXM (K-NOR ) , C ) -EB(IA)* ENTCN 
1 ( RXM (J-2) , BE(IA) ,RXM(K-NOR) , 0) + ( ( EB ( I A ) -E A ( I A ) ) / ( BE ( I A )-A E ( I A ) ) 
2) *(ENT1N(RXM( J-2 ) ,AE( IA) , RXM ( K-NOR ) ,0) - ENT 1 N { RX M { J-2 ) , B E ( I A ) , 

3 RXM { K-NCR ) ,0) ) 

14 AA{ J,K) = AA( J,K ) + AAA 

15 CONTINUE 

155 AA ( J , NOG 1 )«1. 

AA ( J , N0C2 ) = 0. 

- I F ( J-N0RM2) 156, 157, 156 

156 X=RXM ( J-2 ) - EPSX 
CALL SLOP E ( X , AAA) 

GO TO 16 

157 X=RXM ( J-2 ) + EPSX 
CALL SLOPE! X,AAA) 

16 A A( J , NOG 3 ) =- AAA + PSLPCJ-2) 

17 00 2C J=NORM3,NOC2 
DO IE K= 1 , NOR 

AA ( J , K ) =0. 

DO 181 I A= 1 , NO I 

AAA = EA(IA)* DI TON( XR( J -NORM 2 ),AE(IA),XR(K), 1) -EB ( I A) *CITCN 

1 ( XR ( J-N0RM2 ) » BE( I A ) , XR ( K ) , 1 ) + ( ( E B ( I A )-EA ( I A ) ) / ( BE ( I A ) -A E ( I A ) ) ) 

2 *(DIT1N(XR( J-N0RM2), AE( IA) ,XR(K) , 1 ) - D I T 1 N ( XR { J-N0RM2 ) , 8 E ( I A ) , 

3 XR( K ) , 1 ) ) 
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181 A A ( J * K ) = AA(J,K) + AAA 
IF { NOK ) 1 8, 18,920 
920 DO 930 18=1 ,NOK 

AAA =GA (IB)* ENTON{ XR ( J-N0RM2 ) » AG( IB) , XR ( K) , 0) - G8 ( IB ) *E NTON 

1 (XR{ J-N0RM2 ), BG( 18) ,XR( K) ,0) +( (GB( ID)-GA( IB) )/(BG( I B ) — A G ( IB) ) ) 

2 * ( ENT1 N ( XR ( J-N0RM2 ) , AG( IB) ,XR(K),0> -ENT1N t XR( J-N0RM2) , BG( IB) , 

3 XR( K) , 0) ) 

930 AA( J ,K ) =AA { J, K ) -AAA 

18 CONTINUE 

I F (NOO-NOR 1 ) 195,185, 185 
185 DO 19 K=NOR 1 ,NOQ 
AA ( J , K ) = 0. 

DO 19 I A= 1 , NOI 

AAA = EA ( I A )* DITONIXRC J-N0RM2) , AE ( I A) , RXM( K-NOR ) , 0 ) - E B { I A ) * 

1 DITONt XR { J-N0RM2 ) , BE( I A ) , RXM { K-NCR ) ,0 ) + ( ( E B C I A ) -E A ( I A ) ) / { BE ( I A 

2 )- AE ( I A ) ) ) * ( DI T 1 N ( XR ( J-N0RM2 ) , AE ( I A) ,RXM( K-NOR) ,0) - DITIN 

3 (XR ( J- NOR M2 ) , BE ( I A) , RXM ( K-NOR ), 0 ) ) 

19 AA ( J, K ) = A A ( J , K ) + AAA 
195 A A { J , N0Q1 ) =XR( J-N0RM2) 

AAC J , N0Q2 ) = 1. 

CALL DEFLECT { XR C J-N0RM2) , AAA) 

20 AA(J,N0Q3) =-AAA+ PDEFT ( J-N0RM2 ) 

- CALL GAUSS2 (N0Q2, .00000001 ,AA, XA,K1 ) 

GO TO (22,2 1 ) ,K 1 

21 WRITE OUTPUT TAPE 4,77 

77 FORMAT ( 16H MATRIX SINGULAR) 

STOP 

22 DO 23 K= 1 , NOR 

23 RY ( K) = X A( K ) 

I F ( NORM ) 26,26,24 

24 DO 25 K=NO Rl ,NOQ 

25 RAMZ ( K-NOR ) = XA(K) 

26 Cl = XAtNOQl ) 

C2 = XA ( N0Q2) 

X=0. 

NOPE =N0P-1 
FNOP = NOP- 1 
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DO 27 J=1,NCPE 
CALL LOAD (X,W(J ) ) 

CALL SHEAR (X,V( J ) ) 

CALL MOMENT ( X, AM( J ) ) 

XI =X+EPSX 

CALL SLOPE ( X 1 , Y 1 ( J) ) 

Y 1 { J ) = Y 1 ( J ) +C1 
CALL DEFLECT (X, Y( J) ) 

Y(J)=Y(J)+C1*X+C2 
XX ( J ) = X 
27 X = C/FNOP + X 
XX (NOP )=X 
' X=X-3.*EPSX 

CALL LOAD ( X , W ( NOP ) ) 

CALL SHEAR(X, V(NOP) ) 

CALL MOMENT (X»AM( NOP) ) 

CALL SLC° E ( X , Y 1 (NOP) ) 

Y 1 ( NO P ) =Y 1 ( NOP ) +C 1 
CALL CEFL ECT ( X, Y( NOP) ) 

Y(NOP)=Y(NOP)+Cl*X +C2 
WRITE CUTPUT TAPE 4,76 

76 FORMAT ( 3X , 1 H J ,7 X , 4HL0AD, 1 OX, 5HSHEAR , 1 OX , 6HMCMENT , 8X , 5HSLCP E , 

- 1 9X, 1CHDEFLECTION ,8X, 1 HX ) 

WRITE OUTPUT TAPE 4 , 7 8 , ( J , W ( J ) , V ( J ) , AM ( J ) , Y 1 ( J ) , Y( J ) , XX ( J ) , J = 1 , 
1 NOP ) 

78 FORMAT! 14, 6E1 5.8) 

IF (NCR ) 2 82 , 282,280 

280 PRINT 281 , (M,RY(M ) , M= 1 ,NOR) 

281 FORMAT ( 3HRY ( , 12 , 2H ) = , E20. 8 ) 

282 IF(N0RM)288,288,283 

283 PRINT 284, (KA,RAMZ(KA ),KA = 1 , NORM) 

284 FORMAT (5HRAMZ ( , I 2, 2H ) = , £20. 8) 

288 PR INT 286, C 1 

286 FORMAT (5X,3HC1*,E20.8) 

PRINT 267, C2 

287 FORMAT (5X,3HC2=,E20.8) 

285 CONTINUE 
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STOP 

END 
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SUBROUTINE lOAO(X,W) 

DIMENSION A (25) ,B (25) ,QA{ 25) ,CB(25) ,XFC( 100) ,FY( 100) , XR( 10) , RY( 1C) 
1 »XM( 50 f AMZ(50) ,RXM( 10) ,RAMZ( 10),AE(25),RE(25),EA(25),EB(25) 
2,GA(25) ,GB(25),AG( 25),BG(25) 

COMMON A, B,QA,QB, XFO.FY, XR ,RY ,XM, AMZ, RXM, R AMZ , AE , B E » E A , E B t NOL , 

1 NOR» NOF, NOM ,NORM » NCI » NOK » GA ♦ GB » AG , BG 
EPS X= 10 .** (-8 ) 

w=o . . 

1 00 5 I- 1 1 NOL 

2 IF(X+EPSX-B(I))3, 3,5 

3 IF( X+EPSX-A( I ) )5» 4,Ij 

4 WW=QA ( I ) + ( ( QB( I )-QA ( I ) )*(X-A(I ) ))/(B(I)-A(I)) 

w=w+ww 

5 CONTINUE 
RETURN 
END 
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SUBROUTINE SHE AR( X , V ) 

DIMENSION A (25 ) ,B (25) ,QA(25) ,QB(25) ,XFC( ICC) , FY (100) , XR ( 1C) , RY (1C ) 
I*XM(5C) , AMZ ( 50 ) ,R XM ( 1 C ) » RAMZ ( 10 ) , AE( 2 5 ) , BE (25 ) , E A( 25 ) , EB ( 2 5 ) 
2,GA(25),GB(25),AG(25)»BG(25) 

COMMON A, B,QA,QB, XFO, FY , XR, RY , XM , AMZ , RXM, R AMZ , AE ,B E , E A , E B , NOl , 

1 NOR , NOF , NOM , NORM ,NCI,NOK,GA,GB,AG»BG 
EPSX=10.**(-8) 

V=0 . 

IF ( NOL ) 7, 7 , 1 

1 DO 6 1=1, NOL 

IF ( X-A( I ) ) 6,2,2 

2 I F ( X+EPSX-B( I ) )4, 4,3 

3 VV= QA( I ) * ( X- A ( I ) ) + ( (GB( I )-QA( I ) ) * ( ( X- A ( I ) ) ** 2 ) ) /( 2 . « ( B ( I)-A ( I ) ) ) 
l-QBl I ) » ( X-B ( I ) ) - ( (Q8( I )-QA ( I ) ) *( ( X-B( I ) ) **2) )/ (2. *(B( I >-A( I ) ) ) 

GO TO 5 

l4 VV=+Q A ( I ) *(X-A( I) )+( (GB( I )-QA(I JO* ( (X-A( I))»*2))/(2.*(B(I)-A(I))) 

5 V=V+VV 

6 CONTINUE 

7 IF(N0F)10, 10, 75 
75 DO 9 N= 1 , NOF 

IF ( X+EPSX— XFO ( N) ) 9,9,8 

8 VV=FY(N) 

v=v+vv ' 

9 CONTINUE 

10 IF (NOR 113,13,105 
105 DO 12 M= 1 ,NOR 

I F ( X+ EPSX-XR ( M ) 112,12,11 

11 VV= R V ( M ) 

V=V+VV 

12 CONTINUE 

13 CONTINUE 

RETURN 

END 
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SUBROUTINE FOMENT ( X, AM ) 

DIMENSION A(25),B( 25) ,QA(25) ,QB(25) ,XFG( 100) ,FY(100) ,XR( 1 C), RY(IO) 
1 « XM ( 5C ) , AMZ ( 50 } »R XM ( 1 0 ) » R AMZ ( 10),AE(25),BE(25),EA(25),EB(25) 
2,GA(25),GB( 25) , AG (25) »BG(25) 

COMMON A, B,QA,QB, XFO , FY , XR , RY , XM, AMZ, RXM , R AMZ , AE , B E , EA , EB , MOL , 

1 NOR, NGF , NOM ,NORM , NOI » NOK , GA, GB» AG » 8G 
EPSX=10.«*(-8 ) 

AM = 0. 

IF(NOL)8,8,2 

2 00 7 1=1 , NOL 

IFIX-Al I ) )7, 3 , 3 

3 IF ( X+EPSX-B(I ) )5, 5,4 

4 AMM = ( +CA(I)/2. )*( ( X-A( I) ) **2 ) + ( (QB( I ) -QA C I ) ) * ( (X-A ( I )> **3 ))/ ( 6 .« ( 
1BII )-A(I) ) )-(OB(I )/2. )*( ( X-B( I ) 1**2 >- (QB( I >-CA ( I ) ) * ( ( X-B (I) )**3> 
2/ ( 6. * ( B ( I ) —A ( I) )) 

GO TO 6 

5 AMM=( + Q A ( I )/2. )*( ( X — A ( I) )**2)+( ( QB ( I ) -QA ( I } }* ( (X-A ( I ) )** 3 ) )/ (6. * 

1 (B( I ) — A ( I ) ) ) 

6 AM= AM+AMM 

7 CONTINUE 

8 IF(NOF) 12, 12,9 

9 DO 11 N=1 , NOF 

- IFC X+EPSX-XFOIN)) 11,11,10 

1 

10 AMM= F Y ( N ) * ( X-XFO (N) ) 

AM = AM + AMM 

11 CONTINUE 

12 IF(NOR) 16, 16,13 

13 DO 15 M= 1 , NOR 
IF(X+EPSX-XR(M) >15,15,14 

14 AMM= R Y ( M ) * ( X— XR( M) ) 

AM=AM+AMM 

15 CONTINUE 

16 IF (NOM 128 ,28, 17 

17 DO 19 K= 1 , NOM 

I F ( X + EPSX-XM( K ) 119,19,18 

18 AMM= AMZ ( K 1 
AM= AM+AMM 
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19 CONTINUE 

28 IFtNORV) 32,32,29 

29 CO 31 KA= 1 ,NORM 

I F C X+EPSX-RXMIKA) )31,31,30 

30 AMM = R AMZ ( KA ) 

AM= AM+AMN 

31 CONTINUE 

32 CONTINUE 
RETURN 
ENO 
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SUBROUTINE SLOPE! X,Y1 ) 

DIMENSION A(25 ) ,8 (25) ,QA(25) , QB C 2 5 ) ,XFC( 100) ,FY(100) , XR( 10) ,RY( 10 ) 

1 , XM(5C) ,AMZ(50) ,RXM( 1 0 ) , RAMZ ( 1 0 ) , AE (25 ) , B E ( 25 ) , E A( 25 ) , EB { 2 5 ) 

2, GA( 25) ,GB(25), AS (25) ,BG(25) 

COMMON A, B,QA, QB, XFO, FY, XR , RY, XM, AMZ , RXM , RAMZ , AE ,B E, EA* EB» NOL , 

1 NOR,NOF,NON, NORM » NO I , NOK , GA , GB , AG, BG 
Y 1 = 0 . 

00 17 I A= 1 , NO I 
IF(X — AE( IA) ) 17, 17 ,2 

2 IF(N0L)5,5,3 

3 00 4 1=1, NOL 

20 EYY = EA( IA)M ( QA ( I ) /2. ) *ENTON( X, AE( I A) ,A( I), 2) + ( { QB ( I ) -Q A ( I ) ) / 

1 (6. * ( e ( I )-A( I ) ) ) ) MENTON (X,AE( I A) , A( I ) , 3 J-ENTON ( X, AE ( I A ) ,B ( I ) , 3 \ ) 

2 -{QB (I) /2. )*ENTON( X, AE( I A) ,B( I) ,2 ) ) +( ( EB ( I A)-EA (IA))/(BE(IA)- 

3 AE ( I A ) ) ) * ( (0A(I)/2.)»ENT1N(X,AE(IA),A(I),2) M ( GB ( I ) -CA ( I ) ) / 

4 ( 6 . * ( E ( I ) — A ( I ) ) ))* (ENT1N( X, AE( I A) , A ( I ) ,3 )-ENT INC X, AE ( I A ) , B ( I ) ,3) ) 

5 -(QB( I )/2. )*ENT1N(X,AE( IA),B(I) ,2) ) 

Y1=Y1 +EYY 

21 EYY= -EB ( IA)*((QA(I)/2.)*ENT0N(X,BE(IA),A(I),2) + ( ( GB ( I ) -C A ( I ) ) / 
1(6.*(B(I)-A(I)))) MENTON( X,BE( I A) , A ( I ) , 3 )- ENT ON ( X, BE ( I A ) , B ( I ) , 3 ) ) 

2 - { QB ( I )/2. ) *ENTON ( X, BE( IA),B(I),2)) -( ( EB ( I A ) — EA ( I A ) ) / { BE { I A )- 

3 AE(IA)))*( (QA(I)/2.)*ENT1N(X,BE( IA),A(I),2) + ( ( CB ( I ) -QA ( I ) ) / 
4'(6.*( B( I )v- A { I ) ) ) ) * ( ENT IN ( X, BE ( I A ) , A ( I ) , 3 J-ENT 1 N ( X, BE ( I A ) , B ( I ) , 3 ) ) 

5 -(Q8(I)/2.)*ENT1N(X,BE( IA),B(I),2)) 

Y 1=Y 1 +EYY 

4 CONTINUE 

5 IF(N0F)8,8,6 

6 DO 7 N= 1 , NOF 

22 EYY = E A { IA)*FY(N) *ENTON(X,AE( IA),XFO(N) , 1 ) -E B ( I A ) *FY ( N ) *ENTON ( X , 

1 BE( I A ) , XFO ( N ) , 1 ) + ( (EB( IA)-EA( IA) )/(BE( IA)-AE(IA) ) )*FY(N)« 
2(ENT1N(X, AE( I A) ,XFO(N ), 1 ) - ENT 1 N ( X, BE ( I A ) , XFO (N ) , 1 ) ) 

Y1 =Y1 +EYY 

7 CONTINUE 

8 IF(NOR) 1 1 , 11, 9 

9 00 10 M= 1 , NOR 

23 E YY= EA(IA)*RY(M)*ENT0N(X,AE(IA),XR(M),1 ) -EB ( I A ) *R Y ( M ) * ENTON ( X , 

1 BE ( IA),XR(M),1) + ( (EB(IA)-EA( IA) )/(BE(IA)-AE( IA) ))*RY(M) * 
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2(ENT1N(X,AE( I A ) ,X R ( M ) , 1 ) - ENT1N(X,BE( I A ) , XR ( M ) , 1 ) ) 

Yl =Y1 + EYY 

10 CONTINUE 

11 I F ( NOR ) 14,14,12 

12 00 13 K* 1 ,NOM 

» 

24 EYY= E A < I A) *AMZ(K )*ENTON( X, AE ( I A) ,XM( K) , 0) -EB< I A) *AMZ (K ) » ENTON ( X , 

1 BE ( I A ) ,XM(*),0) + ( (EB(IA)-EA(IA) ) /(BE( I A )-*AE ( I A ) ) )• AMZ(K) * 

2 (ENT1N(X,AE( IA), XM(K) ,0)- ENT1 N t X* BE ( I A ) , XM ( K ) , 0) ) 

Yl=Yl+EYY 

13 CONTINUE 

14 I F ( NORM ) 17,17,15 

15 DO 16 KA= 1 , NORM 

25 EYY=E A ( IA ) *RAMZ <K A ) *ENTON ( X , AE ( I A ) , RXM (KA ) , 0 ) -EB ( I A ) * RAMZ ( K A ) * 

1 ENT0N(X,BE(IA),RXM(KA),0) + < < EB ( I A) -E A( I A ) ) / < BE( I A ) -AE( I A ) ) ) • 
2RAMZIKA)* (ENT IN (X ,AE( IA) ,RXM(KA) ,0)- ENT 1 N ( X, BE ( IA ) , RXM ( KA ) ,01) 

Yl =Y1 + EYY 

16 CONTINUE 

17 CONTINUE 
IF(NCK)595, 595,500 

500 DO 59C 16=1, NOK 

IF (X —AG (IE) 1590,590,509 

m 

509 IF (NCL 1530, 530,5 10 

510 00 52 C 1=1, NOL 

26 EYY = GA( IB)*(QA( I ) *FTON ( X, AG ( I B ) ,A(I ),1 ) ♦ (<QB(I) -Q A ( I ) ) / 

1 ( 2. * ( B ( I ) — A ( I 1 )) )*(FTON(X,AG( IB) ,A(I ) , 2 l-FTON ( X ,AG ( I B ) , B ( I ) ,2)1 

2 -Q8( I )*FTON( X»AG ( IB ) ,B( I ) , 1 ) )• + ( ( GB ( IB )-GA { I B) )/ < BG( IB )-AG ( I B ) ) ) 

3 *(QA( I )*FT1N(X,AG( IB) ,A(I) , 1 ) + ( ( QB ( I ) -Q A (I) ) / ( 2 . * ( B ( I ) - A ( I ) ) ) ) 
4* ( FT1N(X,AG( IB), A ( I ) ,2) -FT 1 N( X, AG ( I B ) , B ( I ) , 2 ) ) - QB ( I ) * FT 1 N 

5 (X,AG(IB ),B(I) ,2) ) 

Yl =Y 1 -EYY 

27 EYY =-GB( IB)*(QA( I ) *FTON ( X , BG ( I B ) , A( I ),1 ) + <(QB(I) -QA(I)) / 

1 (2.*(B(I)-A( I))) )*( FT0N(X,BG(IB),A(I),2)-FT0N(X,BG(IB),B(I),2)) 
2-QB ( I )* FTONI X,BG (IB),B(I),1))- ( ( GB ( I B ) -GA ( I B ) ) /( BG( IB)-AG( IB) ) ) 
3* (QA ( I ) * FT1N(X,BG(IB),A(I)»1 )+ ( ( QB( I )-QA ( I ) )/ (2. * (B ( I )-A ( I ) ) ) ) 
4*( FT1N(X,BG( IB), A(I ) ,2)- FT 1 N( X, BGII B) , B ( I ) , 2) ) - QB( I )* FT1N 
5 (X,BG( ie ), B( I ) ,2 ) ) 

Yl =Y1 -EYY 
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520 CONTINUE 

530 IF(NOF)56C»560»540 

540 DO 55 C N=1,N0F 

28 EYY=GA ( IB ) *FY ( N )* FTON ( X , AG ( I B ) , XFO ( N ) , 0 ) -G B ( I B ) * FY ( N ) * FTON(X, 
1BG( 18) ,XFO(N) ,0) + ( (GB,( I B > -G A ( 18) ) / ( BG ( IB ) -AG ( I B ) ) ) * FY ( N ) * 
2 ( FT 1 N ( X, AG( IB) ,XFO(N) ,0) - FT1 N ( X, BG ( I B ) ,XFC <N ) , C ) ) 

Y 1 =Y1 -EYY 
550 CONTINUE 
560 IFINOP )59C,590,570 
570 DO 58 C h‘ * 1 » N 0 R 

29 EYY =GA (IB) *RY ( M) * FTON ( X , AG ( I B ) , XR (M ) ,0 ) -GB ( I B ) * RY ( N ) * FTON(X, 
1 BG ( IB)tXR(M).O) + ( (GB(IB)-GA(IB) ) / ( BG ( I B ) -AG ( IB ) ) ) * RY ( M ) * 

2 ( FT1N(X,AG(IB),XR(M),0) - FT 1 N ( X, BG ( I B ) , XR ( M ) , 0) ) 

Yl =Y 1 -EYY 
580 CONTINUE 
590 CONTINUE 
595 CONTINUE 
RETURN 
END 
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• SUBROUTINE DEFLEC T ( X , Y ) 

DIMENSION A(25) ,B (25) ,QA(25) ,QB(25),XFO( 100) ,FY( 100),XR( 10)»RY( 10) 
1,XM(50),AMZ(50) »RXK(10),RAMZ( 10 ) , AE ( 2 5) , BE (25 ) , EA( 25 ) , EB( 25 ) 
2»GA(25)»GB(25)»AG(25),BG(25) 

COMMON A,B»QA,QB, XFO , FY , XR ,RY , XM, AMZ, RXM , R AMZ , AE ,B £ , EA , EB , NOL , 

1 NOR,NOF,NOM,NORM,NOI,NOK,GA,GB,AG,BG 
Y=0. 

DO 17 I A = 1 » NO I 
IF(X-AE(IA) )17, 17,2 

2 I F ( NOL 15,5,3 

3 DO 4 1=1, NCL 

20 EYY = EA( IA)*((QA(I)/2.)*DIT0N(X,AE(IA),A(I),2) +( { QB ( I )-QA ( I ) )/ 

1 (6.*(B( I )-A( I ) ) )) *(DITON(X,AE( I A) ,A( I ),3)-DIT0N( X, AE ( I A > , B ( I ) , 3 ) ) 

2 -(QB(I)/2.)«DIT0N(X,AE(IA),B(I),2)) + ( ( EB ( I A 1 -EA ( I A )) / ( BE ( I A ) - 

3 AE( IA) ) )*{ (QA( I)/2.)«DIT1N(X,AE( IA),A{ I 1,2) + ( (QB { I ) -QA ( I ) ) / 
4(6.*(B(I)-A(I))))*(DIT1N(X,AE(IA),A(I),3)-DIT1N(X, AE(IA), 6(11,3)) 
5-(QB( I )/2. )*DIT1N(X,AE( I A),B( I ) ,2) ) 

Y=Y + E YY 

21 EYY= -EB(IA)*((QA(I)/2.)*DIT0N(X,BE(IA),A(I),2) +( ( QB ( I ) -Q A ( I ) ) / 

1 (6.*( B( I )-A( I ) ) ))*(DlTON(X,BE( I A) , A(I ) , 3 ) -DITON ( X, BE ( I A ) , B ( I ) , 3 ) ) 

2 -(QB ( I )/2. )*DITON(X,BE( IA),B( I ) ,21) -((EB(IA)-EA(IA))/(BE(IA)- 

3AE( I A ) ) )*( (QA( I)/ 2. )*0IT1N(X,BE(.IA) , A ( I ) ,2) + ( (QB( I )-QA( I ) ) / 
4(6.*(E(I)-A(I))))*(DIT1N(X,8E(IA),A(I),3)“DTT1N(X, BE(IA),B(I),3)) 

5 - < QB ( I )/2. )*DIT1N(X,BE( IA),B( I ),2) ) 

Y=Y+EYY 

4 CONTINUE 

5 IF(N0F)8,8,6 

6 DO 7 N= 1 , NOF 

22 EYY= EA( IA)*FY(N) *DITON(X,AE( IA) ,XFO( N), 1 ) -EB ( I A) *FY( N) *D ITON ( X , 

1 BE( IA) ,XFO(N), 1) + ( ( EB ( I A ) — E A ( IA) )/(BE( I A ) — AE ( IA ) 1 )*FY(N)* 
2(DIT1N(X,AE(IA),XF0(N),1 ) - DIT IN ( X ,B E ( I A ) , XFO( N ) , 1 )) 

y=y+eyy 

7 CONTINUE 

8 IF ( NOR 111,11,9 

9 DO 10 M=1 , NOR 

23 EYY= EA( IA)*RY(M) *DITON( X, AE( IA) ,XR(M), 1 )-EB( IA)*RY(M1* DITON(X, 

1 BE ( I A ) , XR ( M ) , 1 ) +((EB(IA) — EA( IA) )/(BE(IA) — AE (IA)))*RY(M) * 
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2(0IT1N(X,AE( I A ) , X R ( M ) , 1 )- DITlN(X,BE(IA),XR(M),l)) 

Y=Y+EVY 

10 CONTINUE 

11 I F ( NON ) 1 4 » 1 4» 12 

12 DO 13 K= 1 , NOM 

24 EYY=EA( IA)*AMZ(K) *D I TON ( X , AE ( I A ) , XM ( K ) , 0 ) -EB( I A ) * AMZ ( K ) * D I TON ( X , 

1 BE ( I A ) t XM ( K ) t 0 ) + ( (EB(I A)-EA(IA) )/(BE( IAJ-AE (IA)) )*AMZ(K)* ! 

2 (DIT1NU, AE( IA), XM{K) ,0 )- DI T 1 N ( X ,BE ( I A ) , XM ( K ) , 0) ) 

Y=Y+EYY 

13 CONTINUE 

14 I F ( NORM ) 17,17,15 

15 DO 16 KA= 1 , NORM 

25 EYY=EA ( I A )*RAMZ(KA ) *01 TON (X,AE{ IA) , RXM ( K A ) ,0 ) - E8 ( I A ) *RAMZ ( KA )* 

1 DIT0N(X,BE(IA),RXM(KA),0) + ( (EB(IA)-EA(IA) )/{BE( I A ) — A E C IA) ) ) * 
2RAMZI KA)* (DITIN (X ,AE( I A) ,RXM( KA),0)- OIT 1 N( X , BE ( IA ) , RXM { KA ) , C ) ) 
Y=Y+EYY 

16 CONTINUE 

17 CONTINUE 
IF(N0K)595,595, 500 

500 DO 59C I B= 1 ,N0K 

IF (X -AG( IB) )590, 590, 509 
509 IF (NCL) 530, 530,5 10 
SIC DO 52 C .1=1, NCL 

26 EYY=GA( IB)*(QA{ I) * ENTON ( X, AG ( IB) , A{ I ) , 1 ) + ( ( QB ( I ) -CA(I) ) / 

1 (2.* (B ( I ) -A ( I ) )) )*( ENTON (X, AG{ IB) ,A( I ) ,2 )-ENTON{X , AG ( IB ) , B ( I ) , 2 ) ) 

2-QB ( I ) *ENTON( X, AG { IB) ,B( I ), 1 ) ) + ( (GB (IB)-GA(IB) )/ CBG( I8)-AG(I8) ) ) 
3*(QA( I)*ENT1N(X,AG( IB), All ),1) + ( ( QB ( I ) -CA ( I ) ) / {2 . * ( B ( I ) -A ( I ) ) ) ) 
4*( ENT IN ( X , AG ( IB ) , A ( I ) , 2 ) -ENT 1 N ( X, AG ( I B ) ,B( I ) ,2 ) ) - GB ( I ) *ENT IN 
5 ( X , AG ( I B ) , B ( I ) ,2 ) ) 

Y=Y-EYY 

27 EYY=-GB( I B ) * ( QA ( I )*ENTCN(X,BG(IB) ,A{ I ), 1 ) + ( ( QB ( I ) -QA ( I ) ) / 

1 (2.*( B( I ) -A ( I ) ) ))*(ENTON(X,BG(IB),A(I ) , 2 )-ENT ON ( X, BG ( I B) , B ( I ) , 2 ) ) 
2-QB (I )* ENTON (X,BG ( IB ) ,B( I ) , 1 ) )- ( (GB( IB)-GA{ IB) )/( 8G(IB)-AG( IB))) 
3*(QA( I)*ENT1N(X,BG( IB ), All), 1 )+ ( (QB( I)-QA(I) )/(2.*(B(I)-A(I )) )) 
4*(ENT1N(X,BG( IB), A C I ) ,2)-ENT1N(X, BG( IB) ,B( I) ,2) ) - QB ( I ) * ENT 1 N 
5 ( X.BG(IB) ,B( I) ,2 ) ) 

Y=Y-EYY 1 
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520 CONTINUE 

530 IF(N0F)560»560»54C 

540 DO 550 N=1,N0F 

28 EYY=G A ( I B ) *FY ( N) * ENTON ( X, AG ( IB) ,XFO(N) ,0) -GB ( IB )*FY( N ) »£NTCN ( X, 
1BG( IB)»XFO(N) ,0) + ( (GB( IB)-GA(IB) ) / < BG ( I B ) -AG ( I 8 )> ) * FY ( N J * 
2(ENT1N(X,AG( IB)tXFO(N) ,0) - ENT 1 N { X,BG ( I B ) , XFO ( N ) , C ) ) 

Y= Y-E YY 
550 CONTINUE 
560 IF(NOP) 590,590,570 
570 DO 580 M=1,N0R 

29 EYY =GA( IB)*RY(M) *ENTCN( X,AG( IB) ,XR(M) ,0) -GB { 1 8 ) * RY ( M ) * ENTONI X, 
1 BG( 18 ), XR ( M ) ,0 ) + ( (GB( IB)-GA(IB) ) / ( BG ( I B ) - AG ( IB ) ) ) * RY(M) * 

2 ( ENT 1N(X,AG( IB).XR(M) ,0) - ENT 1 N( X, BG ( I B ) , XR ( f' ) ,0) ) 

Y=Y-E Y Y 
580 CONTINUE 
590 CONTINUE 
595 CONTINUE 
RETURN 
6N0 
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FUNCTIGN FTON(X,A,B,N) 
FTON = 0. 

I F ( X- A ) 3 » 3» 1 

1 IF(X-E)3»3,2 

2 FTON= (X-B )**N 

3 RETURN 
END 



FUNCTION FT IN (X,A,B,N) 

FT 1 N =0. 

IF ( X-A ) 3 » 3 V 1 

1 IF ( X-B ) 3 , 3 1 2 

2 FT 1 N = ( X— B )** ( N+ 1 ) + (B-A)*( (X-B)**N) 

3 RETURN 
END 
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FUNCTION ENTCN (X,A V B,N) 

FN=N 

ENTON = 0. 

IF(X-A)5,5, 1 

1 IF(X-e>5,5,2 

2 I F ( A-8 ) 3 » 3 » 4 

3 ENTON = ( X-B ) **( N+1)/(FN+1.) 

GO TO 5 

4 ENTON = { X-B )**{N + 1 ) / { F N+ 1 . ) - ( A— 8 ) * * ( N + 1 ) / ( F N ♦ 1 . ) 

5 RETURN 
ENO 



FUNCTION ENTlNt X, A, B,N) 

FN=N 

ENT1N = 0. 

IF(X-A)5,5, 1 

1 IF(X-B)5»5»2 

2 I F ( A-B ) 3 , 3t 4 

3 ENT1N » { X-B ) **( N+2)/(FN + 2. ) + ( B-A )*( X-B )** (N+1) /( FN+1 . > 

GO TO, 5 

4 ENT1N = ( X-B ) ■**( N+2 } / ( FN+2 . ) + ( B-A) *(X-8)**{N+D / CFN+1.) 

I -(A-B)**(N+2)/(FN+2. > -(< B-A )*( A-B )**( N+ 1 ))/( FN + 1 . ) 

5 RETURN 
ENO 
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FUNCTION DITON(X,A,B, N) 

FN = N 

DITON =0. 

IF(X-A)5,5 f 1 

1 IF(X-E)5,5,2 

2 IF(A-E)3t3»4 

3 DITON = (X-B)*»(N+2) /( (FN + 1 . )*( FN+2. ) ) 

GO TO 5 

4 DITON = (X-B)**(N+2)/((FN+l.)*(FN+2.) ) - (X* ( A-9 ) ** ( N+ 1 ) ) / ( FN + 1 . ) 

l-(A-B)**(N+2)/((FN+l.)*(FN+2.)) + ( A* ( A-B) ** ( N+l ) ) / ( FN+ 1 . ) 

5 RETURN 
END 



FUNCTION DIT1N(X ,A,8,N) 

FN=N 

DIT1N=0. 

I F( X-A) 5t 5» 1 

1 I F ( X- E ) 5 » 5 1 2 

2 IF(A-E)3,3 t 4 

3 DITIN =tX-B)**(N+3) /( (FN+2. )*( FN+3. ) ) + ( (B-A)*( X-B) *»(N+2 ) 1 / 

1 ( ( FN+ 1 . ) * C FN+2. ) ) 

GO TO 5 

4 DITIN = (X-B)**(N + 3)/( (FN + 2. )*( FN+3. ) ) + ( ( B-A ) » (X-B 1 ** (N + 2 ) ) / 

1 ( (FN+1. )*( FN+2. ))-(X*-( A-B) *«(N+2> ) / ( FN+2 . )- ( X* ( B-A J * ( A-B ) ** ( N+ 1 ) ) / 

2 (FN+1.) - (A-B)**(N + 3)/( (FN+2. )*(FN+3. ) ) - ( B-A )* ( A-B )** (N+2 ) / 

3 ( (FN + 1 . ) * ( FN+2 . ) ) + (A*(A-B)**(N+2))/(FN+2.) + ( A * ( 8-A ) * ( A-B ) ** 

4 (N+l ))/( FN+1 . ) 

5 RETURN 
END 



I 
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SUBROUTINE G AUSS2 ( N, E P, A , X, KER 1 
DIMENSION A ( 22, 23 1 ,X(22) 

NPM=N+1 

10 DO 34 L= 1 » N 
KP*0 

Z = 0. C 

00 12 K = L » N 

IF(Z-ABSF(A(K,L) 1)11,12,12 

11 Z-ABSF ( A { K, L ) ) 

K P=K 

12 CONTINUE 
IF(L-KP) 13,20,20 

13 DO 1H J=L,NPM 
Z=A(l, J) 

A ( L, J ) ® A ( KP , J 1 

14 A ( KP , J ) =Z 

20 IF(AESFIA(L,L) 1-EP 150,50, 30 

30 I F { L-N 131,40,40 ' 

31 LP 1 =L+ 1 

DO 34 K = LP 1 , N 
IF(A(K,L) 132,34, 32 

32 RATIO=A(K,L)/A(L»L) 

DO 32 J= LP 1 » NPM 

33 A(K, J)=A(K, J)-RATIO*A(L, J) 

34 CONTINUE 

40 DO 42 1=1, N 
1 1 =N+ 1-1 

DO 43 J=l,1 
JPN= J+N 
S=0.0 

I F ( I I-N141 ,43,43 

41 IIP1=II+1 

DO 42 K=IIP1,N 

42 S=S+A( II ,K J*X(K) 

43 X ( 1 1 ) =(A(II,JPNJ-S)/A(II,II) 

KER S 1 

75 RETURN 



6 ? 



50 K E R=2 
RETURN . 
ENO 
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